A fast direct solver for the advection-diffusion equation using low-rank approximation

Today’s most powerful supercomputers make use of many thousands or millions of parallel processors to achieve high speeds. However, some of the core solution algorithms in computational physics codes are not well suited to this environment. Linear scaling is generally lost as the number of processors becomes very large. New solution algorithms specifically designed for massively parallel execution are urgently needed. In this context, hierarchical algorithms such as the Fast Multipole Method (FMM) [1], H matrices [2] and H matrices [3] have considerable potential. They can solve elliptic PDEs with N degrees of freedom in O(N logN) or O(N) operations and demonstrate excellent parallel scaling [4, 5]. There is one catch, however: hierarchical methods are only useful for PDEs with a Green’s function that decays with increasing distance between two points. Then one can compute a low-rank approximation of the inverse of the system matrix by hierarchically compressing blocks far from the diagonal. Hyperbolic PDEs do not satisfy this criterion. However, Bull et al. [6] showed that by discretising the 1D periodic advection-diffusion equation in time such that the advective term was explicit and the diffusion term implicit, then the linear system was suitable for low-rank approximation. Simple non-hierarchical low-rank approximation methods were able to obtain linear scaling of CPU time with N on a serial processor. In this talk, the ideas are extended to two dimensions and non-periodic domains. New analysis and results are presented. Figure 1 shows the exact solution (left) and solution obtain with the direct solver (right) for propagation of a Gaussian initial condition on a non-periodic domain.